  function oddball(inDir,outDir)

% estrarre i trial per ogni soggetto
% per ogni trial calcolare PLI:
%    -PLI mediando su tutte le coppie di canali
%    -PLI mediando su coppie di canali appartenenti ad una zona (intra)
%    -PLI mediando su coppie di canali appartenenti a zone diverse (inter)
% per ogni soggetto calcolare le medie di PLI sui trial

%inDir='C:\Users\matthew\Desktop\Projects\Oddball\PhaseAnalysis\oddballData
%outDir='C:\Users\matthew\Desktop\Projects\Batch\analysis\oddball\';

pliXBand(inDir,outDir)
pliStat(outDir,outDir);

function pliXBand(inDir,outDir)
bands={[4 8],[8 13],[13 22],[22 49]};
bandLabels={'theta','alpha','beta','gamma'};
nBands=length(bandLabels);

filter=fullfile(inDir,'*.ds');
files=dir(filter);
nSubjs=length(files);

for i=1:nSubjs
    subjName=files(i).name;
    trials=trial(fullfile(inDir,subjName))
    for j=1:nBands   
        
        bf_trials=selectBand(trials,bands{j})
        pliXT=pliXTrial(bf_trials)
        
        outfile=strcat('pli_',bandLabels{j});
        outfile=strcat(outfile,'_');
        subjName=strtok(subjName,'.');
        outfile=strcat(outfile,subjName);
        
        save(fullfile(outDir,outfile),'pliXT');
    end    
end

function [trials]=trial(filename)

% Band filtering
cfg= [];
% Set data
cfg.dataset= filename;
% Remove bad channels, set the trigger, set pre and post-stimulus interval
cfg.channel = {'MEG','-MLT22','-MLO42','-MLO43','-MRT22','-MRO42','-MRO43'};
cfg.trialdef.eventtype= 'OddAudio';  
%cfg.trialdef.eventvalue= eventValue;
cfg.trialdef.prestim= 1.5;
cfg.trialdef.poststim= 1.5;
%Call definetrial to compute trials

trials = ft_definetrial(cfg);
trials.continuous='yes';
trials=ft_preprocessing(trials);
function [bf_trials]=selectBand(trials,band)

cfg = [];
cfg.channel= {'MEG','-MLT22','-MRT22'};
cfg.bpfilter='yes';
cfg.bpfreq=band;
bf_trials= ft_preprocessing(cfg,trials);
function [pli3d]=pliXTrial(bf_trials)

nTrials=length(bf_trials.trial);
nChs=size(bf_trials.trial{1},1);
pli3d=zeros(nTrials,nChs,nChs);
for i=1:nTrials
   pli3d(i,:,:)=pli(bf_trials.trial{i}');  
end


function pliStat(inDir,outDir)

mlc=1:15;
mlf=16:31;
mlo=32:41;
mlp=42:50;
mlt=51:71;

mrc=72:86;
mrf=87:102;
mro=103:112;
mrp=113:121;
mrt=122:142;

%mz=143:151;
labels={'mlc','mlf','mlo','mlp','mlt','mrc','mrf','mro','mrp','mrt'};
areas={mlc;mlf;mlo;mlp;mlt;mrc;mrf;mro;mrp;mrt};

bandlabels={'theta','alpha','beta','gamma'};

for i=1:1%length(bandlabels)

    filter=strcat('pli_',bandlabels{i});
    filter=strcat(filter,'*01*')
    subjs=dir(fullfile(inDir,filter));
    for j=1:length(subjs)
        
        subjName=subjs(j).name;
        load(fullfile(inDir,subjName),'pliXT');
        one_pliXT=pliXT;
        
        subjName(end-4)='2';
        load(fullfile(inDir,subjName),'pliXT');
        
        two_pliXT=pliXT;

        statsingle=statOverTr(one_pliXT,two_pliXT,areas,0.05);
        outfile=subjName(1:end-6);
        [s,outfile]=strtok(outfile,'_');
        outfile=strcat('stat',outfile);
        save(fullfile(outDir,outfile),'statsingle');
    end
     statOverSubjs(inDir,outDir,bandlabels{i},0.01);
end







function PLI=pli(a)

% a is a filtered multichannel signal (time x channels)
% hilbert(a) calculates analytic signal (complex valued) of each
% column of a. Phase Lag Index between channel i and j averaged over
% time bins is stored in PLI(i,j)
% number of channels
N=size(a,2);
PLI(1:N,1:N)=0;
complex_a=hilbert(a);
for i=1:N-1
   for j=i+1:N
        PLI(i,j)=abs(mean(sign(imag(complex_a(:,i)./complex_a(:,j)))));
   end
end
PLI=PLI+PLI';



function [res]=statOverTr(one_pliXT,two_pliXT,areas,alphalevel)

nTrial=min(size(one_pliXT,1),size(two_pliXT,1));
one_trialXintra_area=zeros(nTrial,length(areas));
two_trialXintra_area=zeros(nTrial,length(areas));

interAreaLen=((length(areas)*length(areas))-length(areas))/2;

one_trialXinter_area=zeros(nTrial,interAreaLen);
two_trialXinter_area=zeros(nTrial,interAreaLen);


s=size(one_pliXT);
indx=ones(s(2),s(3));
indx=indx-tril(indx);
indx=logical(indx);

one_avg_pliXtrial=zeros(nTrial,1);
two_avg_pliXtrial=zeros(nTrial,1);


for i=1:nTrial
   
   tone=squeeze(one_pliXT(i,:,:));
   ttwo= squeeze(two_pliXT(i,:,:));
    
   one_trialXintra_area(i,:)=shortPLI(tone,areas);
   two_trialXintra_area(i,:)=shortPLI(ttwo,areas);
   
   one_trialXinter_area(i,:)=longPLI(tone,areas);
   two_trialXinter_area(i,:)=longPLI(ttwo,areas);
   
   one_avg_pliXtrial(i)=mean(tone(indx));
   two_avg_pliXtrial(i)=mean(ttwo(indx));

end

res.intra={one_trialXintra_area,two_trialXintra_area};
res.inter={one_trialXinter_area,two_trialXinter_area};
res.avg={one_avg_pliXtrial,two_avg_pliXtrial};

[res.h_intra,res.p_intra]=ttest(one_trialXintra_area,two_trialXintra_area,alphalevel);
[res.h_inter,res.p_inter]=ttest(one_trialXinter_area,two_trialXinter_area,alphalevel);
[res.h_avgpli,res.p_avgpli]=ttest(one_avg_pliXtrial,two_avg_pliXtrial,alphalevel);


function statOverSubjs(inDir,outDir,bandName,alphalevel)

filter=strcat('stat_*',bandName);
filter=strcat(filter,'*');

files=dir(fullfile(inDir,filter));
nSubjs=length(files);

load(fullfile(inDir,files(1).name),'statsingle');

nAreas=size(statsingle.intra{1},2)
nPairAreas=size(statsingle.inter{1},2);

one_SubjXintraA=zeros(nSubjs,nAreas);
two_SubjXintraA=zeros(nSubjs,nAreas);
one_SubjXinterA=zeros(nSubjs,nPairAreas);
two_SubjXinterA=zeros(nSubjs,nPairAreas);
one_SubjAVG=zeros(nSubjs,1);
two_SubjAVG=zeros(nSubjs,1);

for i=1:nSubjs
    load(fullfile(inDir,files(i).name),'statsingle');

    one_SubjXintraA(i,:)=mean(statsingle.intra{1},1);
    two_SubjXintraA(i,:)=mean(statsingle.intra{2},1);
    one_SubjXinterA(i,:)=mean(statsingle.inter{1},1);
    two_SubjXinterA(i,:)=mean(statsingle.inter{2},1);
    one_SubjAVG(i)=mean(statsingle.avg{1},1);
    two_SubjAVG(i)=mean(statsingle.avg{2},1);
    
    [r,subjName]=strtok(files(i).name,'-');
    subjName=subjName(2:end);
    subjName=strtok(subjName,'.');
    multi.subjNames{i}=subjName;  
end

multi.intra={one_SubjXintraA,two_SubjXintraA};
multi.inter={one_SubjXinterA,two_SubjXinterA};
multi.avg={one_SubjAVG,two_SubjAVG};

[multi.h_intra,multi.p_intra]=ttest(one_SubjXintraA,two_SubjXintraA,alphalevel);
[multi.h_inter,multi.p_inter]=ttest(one_SubjXinterA,two_SubjXinterA,alphalevel);
[multi.h_avgPLI,multi.p_avgPLI]=ttest(one_SubjAVG,two_SubjAVG,alphalevel);

outfile=fullfile(outDir,strcat('statMultiSubjs_',bandName));

save(outfile,'multi');





function [v]=shortPLI(pli,areas)

v=zeros(1,length(areas));
for i=1:length(areas)
    m=pli(areas{i},areas{i}); 
    indx=ones(size(v));
    indx=indx-tril(indx);
    indx=logical(indx);    
    v(i)=mean(m(indx));
end

function [interArea]=longPLI(pli,areas)

v=zeros(length(areas),length(areas));
for i=1:length(areas)-1
    for j=i+1:length(areas)
        m=pli(areas{i},areas{j});
        s=size(m);
        v(i,j)=mean(reshape(m,s(1)*s(2),1));
    end
end

indx=ones(size(v));
indx=indx-tril(indx);
indx=logical(indx);
%indx=logical(indx);%%warning if phase is 0
interArea=v(indx)';


